Surface Melting of the Vortex Lattice in Layered Superconductors: 

Density Functional Theory 
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We study the effects of an afc-surface on the vortex-solid to vortex-liquid transition in layered 
superconductors in the limit of vanishing inter-layer Josephson coupling. We derive the interaction 
between pancake vortices in a semi-infinite sample and adapt the density functional theory of freezing 
to this system. We obtain an effective one-component order-parameter theory which can be used to 
describe the effects of the surface on vortex-lattice melting. Due to the absence of protecting layers 
in the neighbourhood of the surface, the vortex lattice formed near the surface is more susceptible to 
thermal fluctuations. Depending on the value of the magnetic field, we predict either a continuous or 
a discontinuous surface melting transition. For intermediate values of the magnetic field, the surface 
melts continuously, assisting the formation of the liquid phase and suppressing hysteresis above the 
melting transition, a prediction consistent with experimental results. For very low and very high 
magnetic fields, the surface melts discontinuously. The two different surface melting scenarios are 
separated by two surface multicritical points, which we locate on the melting line. 



I. INTRODUCTION 

The melting of the vortex lattice is one of the most 
remarkable aspects of the phenomenology of the high- 
temperature superconductor o 1 ! 2 ! 3 ! 4 ! 5 ! 6 . Such a melting 
transition, in common with other discontinuous phase 
transitions, should involve metastable phases such as an 
undercooled liquid and an overheated solid. Experiments 
on the mixed phase in the layered high-T c superconduc- 
tor BSCCO however indicate an asymmetric hysteresis 
at the melting transition^, involving a supercooled liq- 
uid phase but no overheated solid. In a recent letted, 
we studied the impact of an a6-surface on the vortex lat- 
tice melting in layered superconductor and we found that 
the surface can have a profound effect on the hysteretic 
behaviour. Our approach combined ideas of density func- 
tional theory with a mean-field substrate model proposed 
in Ref. |tj The present paper provides detailed deriva- 
tions of the results outlined in Rcf. 8. In addition, wc 
present several new results, such as a detailed analysis of 
the solid-liquid interface and of the multi-critical point 
at low magnetic fields. 

The absence of a metastable solid phase in ordi- 
nary solids has been the object of great interest in the 
literatureifl*ii'i2ii&±iii&. It is now widely accepted that 
such asymmetric hysteresis can be related to the action 
of surface (pre-)melting. Surfaces act as a natural nuclc- 
ation point for the liquid phase and inhibit the appear- 
ance of the solid at temperatures above melting. This 
is due to the fact that the atoms at the surface experi- 
ence a weaker stabilizing potential. Indeed, the stability 
analysis for a semi-infinite atom chainii shows that the 
enhanced motion of the first atoms makes the 'surface' 
unstable at a lower temperature with respect to the bulk. 
This suggests that the surface may represent a favorable 
nucleation site for the liquid phase. However, criteria 
based on the stability of the solid phase cannot address 



the state of the system beyond the surface instability 
point. 

A comprehensive description of the melting transition 
in the presence of a surface requires the inclusion of the 
effects of both the solid and the liquid phases on an equal 



footi 



.14,16,17 



On a phenomenological level, this can 



be done within the framework of Landau theory using 
the order parameters of the melting transition, i.e., the 
Fourier components of the particle density. The destabi- 
lizing effect of the surface can be accounted foniSii^ by a 
surface term favoring the appearance of the liquid phase. 
Following the discussion in Ref.0, as the simplest exam- 
ple of such a theory one considers a Landau free energy 
of the type 



POO 

F[fi\ = dz 
Jo 



|(^)" + /(m)+^)/iM], (1) 



where z is the distance from the surface, /i the (z- 
dependent) order parameter and f{\i) and the bulk 

free energy and surface destabilizing potential respec- 
tively. For a bulk first-order phase transition with an 
unstable surface we can choose a form 



f(ji) = (a/2)^ 2 - (V3)^ 3 
A(/x) = ( ai /2) M 2 , 



where the coefficients ai,6, c are chosen to be strictly 
positive and temperature independent, whilst tempera- 
ture enters in a(T) via the usual Ansatz for a first-order 
phase transition a(T) — a m cx T—T m , with T m the melting 
temperature and a m — 2b 2 /9c. 

The bulk free energy / shows two minima: one at /ij = 
(liquid) with energy f\ = and a second one at /x s > 
(solid), whose energy depends on temperature. The 
liquid is the stable minimum above the bulk transition 
temperature T m , while below T m the solid becomes the 
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FIG. 1: Top: schematic phase diagram of the pancake vortex 
system. The multicritical point (MC) separates the portion 
of the melting line where the surface melts continuously, O2 
transition (left), from the region where the surface melts dis- 
continuously, Oi transition (right). Lower left: configuration 
of the order-parameter profile fi(z) for an O2 transition at the 
melting temperature T m . In the language of Landau theory 
we use in the paper, the problem maps to the determination 
of the minimal energy configuration of an elastic string, /j,(z), 
subject to a three-dimensional potential landscape, which is 
described by the bulk free energy, plus an additional surface 
term, which pushes the tip of the string towards the liquid 
minimum. In an O2 transition the surface promotes the for- 
mation of the liquid phase upon approaching T m from below. 
The order-parameter profile reaches the liquid minimum con- 
tinuously at the surface for T /* T m , while it remains in the 
'solid' well in the bulk. The liquid propagates continuously 
from the surface into the bulk, thus removing the solid phase. 
Lower right: configuration of the order-parameter profile in 
an 0\ transition at melting. The order parameter reaches a 
finite value on the surface at T m . Melting occurs via a finite 
jump for all z and the solid phase can be overheated. 



stable phase. At the transition the two minima assume 
the same energy and the two phases can coexist. 

In the absence of a surface term, the stable config- 
uration is given by a constant order parameter profile 
throughout the sample. However, the surface free en- 
ergy f\ breaks the invariance of the system along z and 
pushes the tip of the order parameter profile towards the 
liquid minimum at \i\ — 0. This term competes with the 
elastic energy term cx (e?/i / 'dz) . A quantitative under- 
standing of the surface destabilization requires to solve 
the saddle-point equation for the free energy (JJJ. The 
Euler-Lagrange equation provides the differential equa- 
tion 



2 <Pji{z) 
dz 2 



df(rtz)) 
d(j,{z) 



(2) 



which has to be supplemented with the boundary condi- 
tion originating from the surface term, 



,dfi(z) 



dz 



2 = 



dfi 



M(z=0) 



a\\i,(z — 0)- (3) 



Two qualitatively different scenarios are possible de- 
pending on the strength of the surface term or, more 
precisely^, on the ratio a 2 /a m . For relatively large per- 
turbations (a 2 /a m ^> 1) the surface assists the forma- 
tion of the liquid at the melting transition. Approaching 
T m from below, the bulk order parameter remains in the 
'solid' well (see Fig. while at the surface the order 
parameter approaches the liquid minimum at the tran- 
sition temperature continuously. The liquid phase then 
propagates from the surface and the overheated solid is 
removed above T m . Thus, even if the order parameter 
in the bulk jumps discontinuously, the surface melts con- 
tinuously; as the transition is second order at the sur- 
face, in Ref. it is referred to as an O2 transition. On 
the other hand, when the surface free energy is weak 
(a 2 /a m <C 1), the order parameter at the surface jumps 
discontinuously to the liquid minimum across the transi- 
tion as in the bulk, although from a reduced value (see 
Fig. ^) . In Ref. 18, this discontinuous surface melting is 
called an 0\ transition since the transition is first order 
at the surface. In this case, the surface does not inhibit 
the appearance of an overheated solid phase. 

In this paper, we determine the effects of the surface 
on the melting of the vortex system, starting from a 'mi- 
croscopic' theory which accounts for the modification of 
the inter-vortex interactions at the surface. Our program 
is to derive an appropriate Landau-type free energy sim- 
ilar to Eq. Q by considering a coarse grained theory 
where the modulation of the thermally averaged vortex 
density plays the role of the order parameter. The Den- 
sity Functional Theory (DFT) of freezing^ provides an 
appropriate tool for such a study as it is based on the 
change of the free energy from a uniform liquid to a mod- 
ulated stato 4 i 17 i 20 i 21 i 22 i 23 . The only input which is needed 
in the analysis is the direct correlator, which we derive 
within the substrate model approach 9 . Restricting our- 
selves to periodic modulations, one obtains a free-energy 
functional in terms of the Fourier component of the pan- 
cake vortex density at the first reciprocal lattice vectors 
of the ordered lattice. Since these are the natural order 
parameters of the melting transition, the DFT approach 
provides a microscopic derivation of a Landau type de- 
scription of the transition similar to Eq. (QJ . 

As a result of our DFT analysis, we find that for 
very low and very high magnetic fields the surface order- 
parameter still undergoes a finite, although reduced, 
jump at the transition (0\). On the other hand, for in- 
termediate values of B, stray magnetic fields destabilize 
the layers close to the surface, leading to a continuous 
transition at the surface (02)- Consequently, for a wide 
range of experimentally accessible fields we expect that 
surface fluctuations are sufficiently strong to assist the 
smooth propagation of the liquid phase into the bulk and 
to prevent the appearance of the metastable overheated 
solid phase; this result is consistent with the experiments 
of Soibel et ali. The surface continuous and discontin- 
uous transitions are separated by a multi-critical point. 
While a precise location of the high-field multi-critical 
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(a) in-plane 

V z ,z (R) 



(b) out-of-plane 

v z ,z(R) 



force 



FIG. 2: Geometry of the model studied in the paper. A 
layered superconductor fills the half space with z > 0. The 
interaction between two pancake vortices is affected by the 
semi-infinite geometry and by the strong anisotropy, leading 
to a strong logarithmic repulsion between pancake vortices 
which reside within the same layer (a) or to a weak logarithmic 
repulsion between pancake vortices in different layers (b). 



point goes beyond the limit of validity of our analysis, we 
have located the low-field multi-critical point by means 
of an analytical solution of the DFT equations. This 
analysis is further confirmed numerically. 

The paper is organized as follows: we begin by cal- 
culating the modified in-plane and out-of-plane pancake 
vortex interactions (Section II) for a semi-infinite system. 
In Section III, we review the classic DFT of freezing and 
derive a specific functional which describes the vortex 
system in a layered superconductor. We briefly describe 
the results of the novel DFT-substrate approach for an 
infinite (bulk) system (Section IV), where other results 
are available for comparison. We turn to the problem 
with the surface in Section V, where we present an ana- 
lytical solution of the problem and show that depending 
on the value of the magnetic field the system exhibits 
either a 'surface non- melting' (Oi) or a 'surface melting' 
(O2) behavior. Finally, in Section I Vf I we confirm the va- 
lidity of our analytical approach by a direct numerical 
solution of the DFT equations. 



II. MODEL 

We consider a semi-infinite geometry with the super- 
conductor filling the upper half-space z > and describe 
the vortex system within the London theory (see Fig. 
0). We model the system as a stack of two-dimensional 
superconducting layers of thickness d s , separated by a 
distance d, and with a penetration depth A s . We work in 
the limit of vanishingly small Joscphson coupling. The 
basic topological defects are pancake vortices 2 ^ 2 ^ which 
consist of vortices with a two dimensional core limited to 
a single superconducting layer. 

The vortex interaction is mediated by currents set up 
by the vortices via the Lorentz force. A finite sheet cur- 
rent j acts on the vortex core, producing a perpendicular 



($ /c)z x j. 



(4) 



To obtain the force between two vortices, one needs the 
sheet supercurrent generated by any of the two at the lo- 
cation of the core of the other one; then, the interaction 
is found by integrating the force (0J on the radial coor- 
dinate back from infinite distance. The circular current 
defining a vortex is induced by the 2ir twist of the phase 
<p of the complex order parameter ip = \il)\e lv , 



cd s /$. 



4ttA 2 



(5) 



The vector potential A screens the action of the driving 
term V92 = n z x R/i? 2 , annihilating its action when an 
entire flux quantum $ is accumulated. This is the case 
for Pearl vortices in a single superconducting film which 
trap an entire flux <!>o at a distance A off = A 2 /d. However, 
in the case of layered systemsSi^SiSL 2 ^, the presence of 
other layers reduces the ability of the vortex to bind a full 
flux quantum and hence the action persists to infinity. 

A central quantity in the discussion of the vortex- 
vortex interaction is the vector potential field A produced 
by a pancake vortex placed in a semi-infinite sample. In 
the London approximation, a pancake vortex placed at z' 
generates a vector potential satisfying the following set 
of equations 



V 2 A- 



1 



V 2 A = 0, 



;Vip S{z - z'), z>0, 

z < 0. (6) 



Note that Eqs. JBJ are symmetric with respect to z <-» z'. 
In © we treat the layered system within a continuum ap- 
proximation, neglecting small modulations in A (of sub- 
leading order ~ d/X) across the layers. Such a description 
is consistent once the layer thickness d s and the mate- 
rial penetration depth A s are replaced by the inter-layer 
spacing d and the bulk penetration depth A 2 = X 2 d/d s , 
respectively. By solving Eq. ©, we obtain the vector 
potential at the point (R, z) inside the superconductor 
produced by a pancake vortex placed at the origin of the 
layer at z', 



A<p(R,z> Q,z' > 0) 



2A 2 



+OO 



dK J^KR) 
2^ K+ 



x{f z - z ,(K)+(3(K)f z+z ,(K)}, (7) 

with f z (K) = cxp(-K + \z\) and (3{K) = (K + - 
K)/(K + +K) (due to the cylindrical symmetry of © , the 
vector potential has only an azymuthal component, i.e., 
along the direction of the unit vector — n z x R/i? 2 ). 
A convenient quantity, which we will use in the following 
discussion, is the total flux $t(z, z') trapped in the layer 
at z due to a pancake vortex placed at a distance z' from 
the surface; Eq. Q provides the result 



d> t( ,, z ') = ^( e Hz - 2 ' l/A + e - (2+z ' )/A ). 



(8) 



4 



In the following, we will make use of these results to 
derive the pancake vortex interaction between co-planar 
vortices (next subsection) and vortices in different layers 
(Sect. Ell). 



A. In-plane interaction 

We first consider two pancake vortices residing on the 
same layer, cf. Fig. E(a). Their interaction includes 
the contributions from both terms in Eq. JSJ), the driv- 
ing source V</? and the vector potential A. While the 
phase produces currents decreasing cx 1/R, magnetic 
screening is effective at scales larger than the penetra- 
tion depth A. To leading order, the vector potential 
can be written in terms of the trapped magnetic flux 
<I>t(z) = &t(z, z) which threads the layer where the vor- 
tex resides: A<f,{R ^> X,z,z) <&t{z)/2irR. Combining 
Eqs. I0J and (|3J), we can extract the in-plane potential 
for two vortices at a distance R by integrating over the 
radial coordinate from £, the coherence length, to R 



SMALL IN-PLANE SEPARATION tf« X 



V Z>Z (R) w -2e Q d ■ 



4 

4- 



i?< A, 



*t(«) 



$ 



R 



In—, A < R, 

A 



(9) 

where &t(z) is the magnetic flux trapped by a pancake 
vortex placed at a distance z away from the surface. From 
(JSJ), we find that the trapped flux $t(z) interpolates be- 
tween the values (d/A)$o on the surface and (d/2A)<i>o in 
the bulk. 



(10) 



However, we can safely ignore the small modifications 
of order d/X in the potential 10 arising from magnetic 
screening and assume an in-plane interaction which is 
independent of z. Within this approximation, vortices in 
the same layer feel a repulsive logarithmic interaction, 



V Z , Z (R) w -2s dln(R/£). 



(11) 



This potential corresponds to the one between charges in 
the two-dimensional Coulomb gas, or, equivalently, the 
one-component plasma (OCP), with charge e 2 — > 2end. 



B. Out-of-plane interaction 

For two pancake vortices in different layers (cf. Fig. 
I2jb)), the current (JjJJ which enters the Lorentz force Q 
does not contain the contribution from the driving phase 
Vyj. Hence, the force in |(5J is due to the vector potential 
A alone. The central quantity in the discussion is the 
vector field J7J which is produced by a single pancake 
vortex placed in a semi-infinite sample. 



(a) semi-infinite 

v z , z -(R) 



(b) bulk 



FIG. 3: For small in-plane separation, i.e., R <C A, the total 
energy for the out-of-plane interaction in a semi-infinite sys- 
tem (a) is equivalent to the one of a translation invariant bulk 
(b) system V*_ Z ,(R). 



Inserting the vector potential J7J into and (0, we 
obtain an expression for the force between two pancake 
vortices residing in different layers z ^= z' . The integra- 
tion over R provides the out-of-plane interaction 



V. 9 ,{R) 



-e d* r d K J ^ 
Jo 



KK+ 

x[f z _ z ,(K)+(3(K)f z+z ,(K)}. 



(12) 



This interaction is given by the sum of a bulk- (first) 
V^_ Z ,{R) and a stray-field term (second) V z s 2 , (R), where 
the latter becomes negligible at a distance ~ A away from 
the surface. 

For small in-plane distances, i.e., R <C A, the contribu- 
tion of the stray-field potential to the overall vortex in- 
teraction energy can be neglected: V ZZ ,(R) <C V Z _ Z ,(R). 
Thus, for R -C A, the potential coincides with the bulk 
one M 



V z>z/ {R)tt So d-{ 



R 2 



4\z- 
R 
A' 



z'\\ 



\z -z'\ < A, 



« « A. 



As expected, the out-of-plane interaction is smaller with 
respect to the in-plane expression i|II|) due to the small 
pre-factor d/X. However, the out-of-plane interaction ex- 
tends over many (X/d) layers, underlying its importance. 
In the regime R -C A, the attractive interaction between 
a single pancake vortex and a semi-infinite stack then is 
reduced close to the surface due to the missing planes for 
z < 0. 

On the other hand, for large in-plane distances, i.e., 
R 3> A, we find relevant modifications of the out-of-plane 
interaction at the surface 

V„,(R)*< So die-"-'<»'ln!- 

+-4- ( * + *'" A ('4 + f)- ("> 
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LARGE IN-PLANE SEPARATION R » X 



(a) semi-infinite 

V---(R) 



(b) bulk 

vt z '(R) 



(c) image vortex 



FIG. 4: For large in-plane separation, i.e., R ^> A, the to- 
tal energy for the out-of-plane interaction in a semi-infinite 
system (a) can be split in (b) a translation invariant bulk 
term V z b _ z ,(R) and (c) an additional one V Z+Z ,(R) that can 
be interpreted in terms of an image vortex placed in —z. The 
stray field contributes via an additional algebraic correction, 
cf. (£3 



The first term has a bulk origin, whereas the terms which 
follow vanish away from the surface. Combining terms, 
we can rewrite (|14f) in the more convenient form 



V z , z ,(R) = 2e Q d 



$ 



A A R 



A 



(15) 



We identify two different contributions: a logarithmic 
attractive interaction (first) and an algebraic repulsive 
potential (second). The logarithmic interaction can be 
analyzed by considering the presence of mirror vortices 
in z < 0. The two terms in <S>t(z, z') can be viewed as 
the contributions of two bulk vortices, one at z and a 
mirror one at — z, see Fig. 0] Hence, if we consider only 
the logarithmic term in (|15|) . the interaction between a 
pancake vortex and a stack in a semi-infinite system is 
equivalent to that in a bulk system. Modifications appear 
due to the stray magnetic field, generating the residual 
algebraic potential in i[T5)l. This correction oc X/R is 
responsible for the weak algebraic interaction between 
vortex tips at large distances^ .. 

Summarizing, for short in-plane distances, R -C A, 
the strength of the out-of-plane potential is not affected 
by the semi-infinite geometry (see Fig. 0J and the effect 
of the surface is limited to the lack of superconducting 
planes at z < 0. On the other hand, for large in-plane 
distances, R 3> A, the missing planes are compensated 
by the contributions from the mirror term in l|14|) . see 
Fig- El Surface effects then are due to the stray magnetic 
field generating an additional algebraic repulsion, cf. Eq. 

03- 



III. DFT-SUBSTRATE APPROACH 

Here, we briefly discuss the DFT-substrate approach 
we will use in our analysis (see Ref. for details). In 



the classical DFT one writes the total grand canonical 
free energy difference from the homogeneous liquid phase 
as a functional of the (averaged) density profile. In an 
anisotropic system such as the pancake vortex system, it 
is convenient to separate the in-plane dependence from 
the out-of-plane one. Hence, the DFT free energy takes 
the form 



sn[ Pz (R)} 



T 

i r dz' 



dz 
a 



p z (R) ln^^-^(R) 



d 2 R<V z (R)c z , z ,(|R-R|)<W(R) , (16) 



where p z (R) is the averaged vortex density in the layer 
z and Sp z (R) = p z (R) — p is the deviation from the 
homogeneous liquid density p. The only input needed 
in the DFT free energy is the direct correlation function 
c z,z'(R)- In our approach we implement the substrate 
model 9 for the pair-correlation function by separating the 
contributions of the strong in-plane logarithmic repulsion 
from the weak out-of-plane but long-range attraction, 



c z , z ,(R) = dc 2D (R)S(z - z') - V Z , Z ,(R)/T, 



(17) 



where V Z . Z >(R) is the out-of-plane interaction of Eq. 112JI : 
the above approximation for the correlator is described 
in detail in Ref. Within the planes, vortices are 

strongly correlated due to the repulsive logarithmic in- 
teractions (we neglect small contributions of order d/X 
and use QJ ). Hence, we can approximate c z>z {R) with 
the direct correlation function c 2D (R) of two-dimensional 
logarithmically interacting particles (also known as one- 
component plasma, OCP). Instead of using an approxi- 
mate scheme, such as the hypernetted chain equation 20 , 
we use results of Monte Carlo simulations of the two- 
dimensional OCP at various coupling constants T = 
2eod/T to extract c 2D (R). On the other hand, the out-of- 
plane direct correlator c z . z >(R) is determined within per- 
turbation theory^ in terms of the weak out-of-plane po- 
tential; neglecting higher-order correlations, we approxi- 
mate it with the leading unperturbed value —V ZjZ '(R)/T. 

At a mean-field level the thermodynamically stable 
state corresponds to the minimal free energy configura- 
tion of the functional (|16(l . Then, the density functions 
p z (R) must obey the saddle-point equation 



hi 



'''- lRl - j^j J d 2 R' c z , z ,(\R-R'\)6 Pz ,(R'). (18) 



A key quantity in our discussion is the molecular 
field 19 ' 32 ! 33 £ Z (R) defined through 



&(R) = ln[ Pz (R)/p]. 



(19) 



Combining the saddle point equation (|18fl with (|19fl . the 
molecular field becomes 

6(R) = J^fJ d 2 R' c ZiZ ,(\R--R'\)5p z ,(R') (20) 

and represents the average potential produced by the 
modulated density. However, whereas i|19|) defines the 
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molecular field everywhere, Eq. <|20|l is only valid at the 
minimum. 

Next, instead of seeking the exact form p z (R) which 
solves the non-linear integral equations l|18f) . we restrict 
our analysis to a simple family of periodic functions which 
describes the modulations of the density in the solid 
phase. In the following, we concentrate on the simplest 
case, retaining only the first Fourier components of the 
density in a triangular lattice, 



J2 /'■■ 

|Ki |=G 



e «i-H = l + ^ ff (R) j (21) 



where the vectors are the first reciprocal lattice vec- 
tors of the frozen structure of length G which is re- 
lated to the lattice constant a A via G — 47r/v / 3aA, and 
|U Z = Sp z (G)/p is the Fourier component of the den- 
sity with wave length G. In the following discussion, 
we consider the vortex system as incompressible and ne- 
glect the relative density change Sp(K = Q)/p, hence 
the length of the first reciprocal lattice vectors can be 
written as G = {8n 2 p/V3) 1/2 (see Ref. E2 for a dis- 
cussion of how to include the effects of a finite com- 
pressibility). The function <?(R) = Y^\Kx\=G e * Kl R = 
2cos(Gx) + 4cos(Gx/2) cos(v / 3G?//2) includes the sum 
over the six first reciprocal lattice vectors. We also write 
a similar Ansatz for the molecular field 



e,(H.) = c+e.s(R) 



(22) 



retaining only the zeroth and first Fourier components, 
£ z and £z respectively. The Fourier components of p(R) 
and £(R) are not independent and can be related through 
(|19fl . Filtering out the zeroth and first Fourier compo- 
nents of Pz(R) = p{l + p z g(R)) = /oexp(( z + £ z g Kl (R)) 
and using 



- [ d 2 Rg(R) = 0, - { d 2 R[g{R)] 2 = 6, 

a J a a J a 

we obtain the dependencies 

Cz = (j-z = *'(6)/6, 

where we have defined the function 



$(£) = In 



- I d 2 Re^ 



(23) 



(24) 



(25) 



Inserting the above Ansatz for the order-parameter pro- 
file Sp into the free energy (|lf>|> . we obtain the associ- 
ated free energy density functional (A is the sample area) 
5ui/T = (1/ pA)5$l/T expressed through p z 



Soj[p z 
T 



where 



D dz 



'°°dz' _ 

66Mz-*(6)-3/ —j-c ZjZ/ (i z ifi z 



c 2D dS(z 



V Z . Z> (G)/T 



(26) 



(27) 



is the dimensionless Fourier transform^ of the direct cor- 
relation function c z z i — c z , z >(G) evaluated at the first re- 
ciprocal lattice vector G. Consistently, we have defined 
c 2D = c 2D (G) . The effects of stray fields are accounted for 
within the out-of-plane interaction V ZtZ >(G) of Eq. (|12|) . 



Vz,z>(G) 
T 



2-KpTd 



,( -G+\z-z'\ , if± 

G 2 \ 2 G+V T G+ 

-a(fz- z < + Pfz+z')- 



G-\- — G 



-G+(z+z') 



+ G 



(28) 



In J2EJ) we have defined f z _— f z (G) = exp(-G+|z|), a — 
a(G) = 2np£d/G 2 X 2 G + , /? = (G+ - G)/(G+ + G), and 
G + = \J G 2 + 1/A 2 . Note the additional factor p due 
to the dimensionless definition of the Fourier transform. 
Combining (|26() with (|27|l and (|28|l and after few simple 
algebraic manipulations we obtain the functional 



Scj\pz 
T 



dz 
1 



+ 



10 

3 a 
6aG 

gTgZ 



dz' 



T 



(fz-z> + fz+z>)(^z - /v) 2 



dz 
~~d 



dz' 
~d~ 



ip-zPz 



(29) 



The first term in (|29(l describes the local two-dimensional 
free energy of a uniform system 



5u>™(p)/T = 6£(p)p - - 3c- M 2 



(30) 



where £ has to be understood as a function of p through 
l)24[l. In 13()|l we have defined the correlator 



K tack (G)/T, 



(31) 



where the out-of-plane interactions contribute to the cor- 
relator through the total stack potential 



Kta^(G) 



-a 



dz 



\fz-z' + fz+z') 



2a 
~dGl 



(32) 



The additional non-local terms in (|29|1 quantify the en- 
ergy cost due to variations of the order parameter (the 
non-local 'elastic' term, second line) and the effect of the 
surface (third line). 

The free-energy (|29|) has to be compared with the phe- 
nomenological Landau-type functional of Eq. Q). Both 
expressions exhibit a similar structure composed by a 
sum of bulk, 'elastic', and surface terms. However, the 
analysis of the free energy i|29[l is complicated by the 
non-local nature of the 'elastic' and surface terms. 



IV. BULK SYSTEM 

Before concentrating on the surface problem, we briefly 
review the case of an infinite bulk system. The results 
in this section will be the basis for the discussion of the 
melting transition in the semi-infinite geometry which we 
will present in the following sections. 
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In a bulk system, we can drop the surface term and 
the free energy (|29[) becomes 



T 



chvSu} 2 ° h (p z 
d 



3a 

T 



T 
dz' 



fz-x>(jJ-z' - HzY 



(33) 



Below, we first concentrate on the thermodynamically 
stable state (see also Ref. l22j) and then analyze the prop- 
erties of a solid-liquid interface. 



A. Thermodynamic phase diagram 

In thermodynamic equilibrium, all superconducting 
layers are equivalent and the order parameter \i z be- 
comes independent of the layer position z, hence we write 
fx z = fi. For a constant order-parameter profile the sec- 
ond term in l|33|) is zero and only the bulk free energy 
of Eq. H30|) is relevant. The key quantity is the effec- 
tive three-dimensional correlator c 2 ° b which is given by 
the sum of two contributions: the OCP term c 2D and the 
stack potential K> tack (G) of Eq. Whereas the first 

depends only on the temperature, the latter depends also 
on the vortex density and thus on the magnetic field, 



C(T,B) =c 2D (T) + 
= c 2D (T) + 



Airpeod 1 
TG 2 1 + A 2 G 2 



(34) 



2ttT [1 + (8ir 2 /V3)B/B x ] ' 



where we use p/G 2 = %/3/(8tt 2 ) and A 2 G 2 = 
(8ir 2 /V3)B/B x , with B\ = $ /A 2 . 

The temperature and magnetic field dependencies en- 
tering c 2 ° b (T,B) change the coefficient of the quadratic 
term as in the </> 4 Landau theory. As a function of p,, 
the free energy exhibits the shape of a Landau theory 
describing a first-order phase transition (this function 
is plotted in Ref. 12^1 : for small values of the correla- 
tor c 2 ° b , corresponding to large temperatures and fields, 
Su 2D (p) exhibits only one minimum at p\ = with a 
value <5w 2D (0)/T = 0, corresponding to a (homogeneous) 
liquid phase. Decreasing the temperature and/or field 
(which corresponds to increasing values of c 2 ° b ), a second 
local minimum at fi s > (metastable solid) appears in 
addition to the liquid one at p\ — 0. Freezing occurs 
when the liquid and solid minima assume the same value 
of the free energy, i.e., when Sui 2D (p s ) = 0. This condi- 
tion is equivalent to a simple equation for the correlator— 

?° h (T,B) = c c w 0.856. (35) 

Solving this equation for B yields the melting line, see 
Ref. Going to even lower temperatures, c 2D fur- 

ther increases, the solid minimum decreases further, 
Sw 2U (n s )/T < 0, and the crystal becomes the only ther- 
modynamically stable phase. 



B. Solid-liquid interface 

Next, we consider a solid-liquid interface in an infi- 
nite system, for which both terms in ll-i-il) are impor- 
tant. Along the melting line B m (T) the solid and the 
liquid assume the same value of the free energy and can 
coexist. Here, we analyze the situation where half of 
the system is in the liquid state and the other half is 
solid. The corresponding profile of p, z takes the form of 
a soliton with boundary conditions p. z ^-oo = (liquid) 
and p z ^+Qo = p s ~ 0.51 (solid), where 5uj 2 ° b (0)/T = 
5uj 2 ^ b (p s ) /T = 0. The properties of the interface be- 
tween the two phases depend on the non-local term in 
the free-energy functional i|33|l . For a large part of the 
phase diagram (see Appendix [XJ, it is possible to ap- 
proximate the full non-local theory with a local one by 
proceeding with a gradient expansioni&Si of the kernel 
in 133fl : inserting p z > w \i z + (dp, z /dz){z' — z) into 13311 . 
we obtain 



T 



00 dz 
.™ d 



P/djiz 
2 V dz 



(36) 



with the clastic scale 



+ OO 



dzf z z 2 = 



3\/3 



12a 

dGl 

r 



2^ [l + (8Tr 2 /V3)B m (T)/B x ] 2 ' 



(37) 



This local approximation describes well the full non-local 
free energy if the profile p z varies slowly over the exten- 
sion 1/G+ of the kernel. We have checked that this con- 
dition is fulfilled for T > 0.04 e <i, which corresponds to 
small and moderate magnetic fields B < Q.5B\, by com- 
paring the shape of the soliton in the non-local and local 
theories (see Appendix [XJ. 

We can easily estimate the main properties of the 
solid-liquid interface, e.g. its width and free energy cost. 
At small magnetic fields, we use the gradient expansion 
and write the Euler-Lagrange equation of the free-energy 
functional (|36|l 



dz 2 



dll; 



(38) 



and the corresponding energy conservation (the integra- 
tion constant is zero for a soliton solution) 



df± = 1 26lj™M 
dz U T 



(39) 



The energy E s \ of the interface (soliton) is given by 

T 



K] = r 1 ^ e ( d ^\ 2 _ r .... 



d V dz 



d Jo 



dp 



T 



(40) 
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FIG. 5: Profile of the order parameter /i 2 describing an in- 
terface between a liquid (^t z ^_oo = fJ-i = 0) and a solid 
(fiz^+oo = /i s ~ 0.51). The shape has been obtained numer- 
ically by solving the differential equation 1381 . The width of 
the interface is approximately £ a ~ 4.05 £ (we define it as the 
full width at half maximum of the function dfi z /dz), where £ 
is defined in 13711 . 



where <W^ x ~ 0.0065 T is the barrier in the energy den- 
sity between the solid and liquid minima on the melting 
line and fj, s ~ 0.51. The constant C is of order unity 
and requires to evaluate the integral in (|40|l : a numerical 
calculation gives C « 0.69. 

Next, we determine the width of the soliton £ c , which 
we define as the full width at half maximum of the 
derivative d\i z jdz. Approximating dji z /dz « jj, s /£ e in 
11401 . we obtain an estimate for the soliton energy E s \ w 
T£ 2 (£ /d)(/j s /£ c ) 2 . Combining this formula together with 
we find 



/'s 



y/25u 2 °jT 



(41) 



From l|41(l . we obtain £ c w 6.33 £; however, a precise de- 
termination of the proportionality factor requires a nu- 
merical study, which yields £ c s» 4.05^?, see Fig. El At low 
fields, B — ► (and T — > T BKT ), the soliton becomes wider 
than the bulk penetration depth, £ c 8A, see l{4"T|) and 
(|37|l . This confirms the validity of the gradient expansion 
for B — > 0, since in this limit 1/G+ ~ A and ^ e G+ w 8. 

On the other hand, for large £?, the gradient expan- 
sion is not valid and we cannot use Eq. I|41[l . However, 
we have determined numerically the shape of the soliton 
(see Appendix lAl using the full non-local theory (|33|l and 
found its width to decrease with increasing B (the full 
theory produces a sharper interface than the one origi- 
nating from the gradient approximation). We can inter- 
pret this results in terms of the intra-layer interactions 
dominating over the inter-layer ones, leading to a system 
of decoupled 2D planes which melt independently. 



superconductor occupies only a half-infinite space z > 0, 
ii) stray magnetic fields modify the vortex interaction 
and hence the direct correlation function c z>z i near the 
surface. 

The configuration of the system is derived from the 
saddle-point equation of (|29|l . Minimization with respect 
to n z provides us with the integral equation 

d 



6a 

to 
12tvG 

G + G 4 



■Ifz 



'dz' - 

d /; 



+ fz+z']{Pz - Hz') 

(42) 



• f'z 



0. 



The surface produces two additional terms when com- 
pared to an infinite system: a mirror term (oc f z + z >) in 
the non-local elastic force and a pure surface term (sec- 
ond line) , cf . the discussion of the out-of-plane interaction 
in Sec. inB 

In the following, we first acquaint ourselves with the 
formalism by studying the B as and large field limits. 
Then, we present the full analysis at low but finite values 
of the magnetic field. 



A. Large magnetic field and B « limits 

Let us consider first large magnetic fields, B ^> B\. 
In this limit, the contribution of the out-of-plane poten- 
tial vanishes (G w (B/^q) 1 ^ 2 — > oo and a — > 0) and the 
system decouples into independent two-dimensional sys- 
tems. Only the intra-layer interaction remains relevant 
in this limit and the free energy is simply given by the 
bulk local contribution (|3U|) . The saddle point equation 



d lta 8 U %Q*,)/T = 



(43) 



does not depend on z and the solution is given by a con- 
stant order parameter profile. Thus, no modifications 
occur at the surface for B ^> B\. 

In the opposite limit of B -> (G -> 0), G+ w 1/A 
and a remain finite. The surface term in l|29() is negligible 
and the saddle point equation 



T 



6a 



°° dz' - - 

—J- [fz-z'+fz+z'\{P-z 



fJ-z 



0, 



contains an additional term as compared to (|43|l . How- 
ever, we can render the problem translational invariant 
on the whole z-axis, by introducing mirror vortices in 
z < via /u_ z = fi z and changing z' — > —z' in the in- 
tegral of the mirror term. As a result, we can map the 
semi-infinite system back to the infinite bulk and the con- 
stant bulk solution remains valid also in this limit. We 
can conclude that for both B sa and B ^ B\ the sur- 
face docs not lead to a modification of the bulk behavior. 



SURFACE MELTING 



B. Low magnetic fields 



The presence of a surface modifies the bulk free energy 
functional l|33|l in two ways as it is clear from l|29|) : i) the 



At finite magnetic fields, all three terms of the saddle- 
point equation (|42|l are relevant. To understand the effect 
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of the surface term (third), it is convenient to look at the 
associated free energy expression (cf. I|29|) L which can be 
written as a weighted square of /i z , 



T 



6aG 
G + G, 



(/; 7 



(44) 



This term favors the appearance of the liquid on the sur- 
face, similarly to the local surface term in 0J. In real 
space, 1)440 is associated with the oc 1/R repulsive poten- 
tial induced by the stray magnetic field, cf. the second 
term in (|15|l . In analogy to the Landau-theory descrip- 
tion of Eq. ||TJ, this DFT functional can yield two dif- 
ferent surface melting scenarios: discontinuous, 0\, and 
continuous, O2 (see introduction). 

To make progress analytically, we have to simplify the 
non-local terms in (|42l) . Concentrating on the bulk, i.e., 
z 3> 1/ G+ , both mirror and surface terms can be ignored. 
For not too large values of the magnetic fields, we can 
adopt a gradient expansion of the non-local elastic term 

fz-z>, cf. Sec. II V Bl and Appendix 1X1 As a result we 
obtain the differential equation 



>d 2 fj, z 



dz 2 djj, z T 



(45) 



with £ 2 = (6a/d) J +o ° dz f z z 2 = !2a/dG 3 +1 cf. Eq. ©. 
Equation l|4*5|) has to be completed with a boundary con- 
dition, which has to be provided by the surface term in 
the integral equation. In the following, we restrict the 
analysis to small values of the order parameter at the 
surface in order to study the boundary between continu- 
ous and discontinuous surface melting scenarios. In this 
regime the bulk potential (|30[l can be approximated as 



6lo™^)/T » 3(1 - c c ) M 2 



(46) 



where we used $(£) w 3£ 2 (cf. (JU), fi = $'(0/6 « 
£ (using l|24l) together with (|23[l ) and the value of the 
critical correlator c 2 ° b — c c « 0.856 at melting. The 
saddle point equation takes the form of a linear integral 
equation 

_ r 00 dz' - 

6(1 - c c )fi z + 6a —[fz-z' + fz+z'](pz - Hz')+ 

\2aG f°°dz' - 
+ ~ : ~ / — /»+»' AV = ( 47 ) 



G + G+Jo d 

with three terms: a potential force (first), a non-local 
elastic force (second) and the drive at the surface (third) . 
Separating the non-local terms from the local ones, we 
obtain 



6(1 - Cc) + 



12a 



dG A 



fj, z = 6a 



00 dz' - 



+ Pfz+z>}Hz> 



(48) 

where f3 was defined in l|28|) . The second term in the left 
hand side carries the normalization of the kernel f z (from 
the integral 6a. J °° (dz' /d)[f z - z i -I- f z + z '] = l2a/G + d in 



the elastic term of Q47|l). Equation (|48|l can be rewritten 

as 



Hz = 



2(1 + r) Jo 



dz'[f z - z , +M*+AV*, ( 49 ) 



with r = dG+(l - c c )/2a = 6(1 - c c )/(G+£) 2 . 

This kind of integral equation is commonly found in the 
study of boundary problems, e.g. in the analysis of sur- 
face effects on the superconducting transitional^. For 
the latter case, the equation which determines the super- 
conducting gap A(r) is a non-local integral equation with 
a structure similar to (f^Hf) (or Close to the super- 

conducting transition T c , a gradient expansion reduces 
the non-local equation into the local Ginzburg-Landau 
equation. The presence of the surface enters the kernel 
of the non-local theory via an additional contribution, 
corresponding to the term oc f z + z i in (|48|l . However, a 
major difference between the two problems is given by 
the different nature of the bulk phase transitions, first 
order for the melting and second order for the supercon- 
ducting transition. In the case of a second-order phase 
transition the coefficient of the quadratic term in the bulk 
free energy goes to zero at the transition point, cf. the co- 
efficient a ~ (T — T c ) in the Ginzburg-Landau equation. 
Hence, the term corresponding to the force oc 6(1 — c c ) 
in the LHS of |@SJ| is absent. In Eq. (03 the difference 
between a first- and a second-order bulk transition enters 
in the normalization of the integral in the RHS, through 
the parameter r, which is r ^ for a first-order transition 
and r = for a second-order one. 

In general, the solution of an integral equation of the 
type l|49l) is a non-trivial task, which usually cannot 
be carried out exactly. However, in the present case a 
straightforward solution is possible due to the particular 
exponential form of the kernel. In fact, taking the second 
derivative of ll^ and making use of the identity 



dz 2 



~ G +\ Z - Z '\ = -2G+S(z - z') + G^e- G +I 2 ~ 2 'l, (50) 



we find that the integral equation yields the same differ- 
ential equation (|45|l as derived previously in the bulk by 
means of a gradient expansion, 



(1 + r) ^S i " £2 S i=r(£G+)2 ^ 

6(1 -c c )m*; (51) 



here, we have used the linearized force d^Suj^dij/T = 
6(1 — c c )/i z and we can drop the small renormalization 
factor (1 + r) w 1 in the LHS of the equation^!. The dif- 
ferential equation (|51fl admits two exponential solutions. 
This has to be compared to the case of a second-order 
bulk transition, where, due to the absence of a linear 
term in the bulk free energy, one obtains the differen- 
tial equation d 2 /i z /dz 2 = which is solved by a linear 
function- 8 . In both cases, the integral equation (|49|l is 
equivalent to the bulk differential equation l|51|l and the 
bulk solution goes through the surface region. Hence, 
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the effect of the surface terms oc j3 is only to provide a 
boundary condition at z — but does not modify the 
bulk solution. This lucky coincidence is due to the par- 
ticular exponential structure of the kernel. Usually, the 
bulk solution is not valid in the vicinity of the surface 
and the problem becomes much more difficult to solve. 

Next, we derive the boundary condition which is pro- 
vided by the surface term. We first write the general 
solutions of the bulk differential equation 

fi z = Ae^+z + Be-~< G + Z , (52) 

where j 2 = r/(l + r); the boundary condition then fol- 
lows from the ratio A/B. Here, we retain the small cor- 
rection (1 + r) in (|51|l as we base our analysis on the 
linearized integral equation (dropping this term leads to 
results which are correct to order r, in agreement with 
the precision of the gradient expansion). Inserting this 
Ansatz back into (|49|l . we obtain the condition 

a(-I 1 -)+bU U=o. 

VI - 7 I + 7/ V1 + 7 1-7/ 
This requirement selects a unique value of the ratio B/A, 
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FIG. 6: Graphical solution of I56H for /io- We plot separately 
the RHS of gSJ, ^/26uj™ b ( f j. )/T, and the LHS, £Gfi , for 
different values of the slope £G. For £G < \J 6(1 — c c ) we find 
an intersection point at a /to > (besides the solution /to = 
0). This finite value is the residual order parameter on the 
surface for the 0\ scenario (surface non-melting) . For £G = 
\J 6(1 — c c ) the straight line is tangent to ^25u!^° b (fio) /T at 
/io = 0, the finite solution has disappeared, and only /io — 
remains (multi-critical point). For £G > \J 6(1 — c c ) the only 
solution is at /io = (O2 scenario, surface melting). 



B _ ^ft+l) + (/?_!) 

A 7 tf + i)-tf-iy 



(53) 



Finally, we can calculate the logarithmic derivative at 
z = 0, which gives the relevant boundary condition for 
our forthcoming analysis, cf. Eq. 0, and we find the 
result (we define \J Z = d z (i z ) 



= G+7 ^L_| = G+ LJL G . (54) 
=0 + ' A + B + l + (3 y ' 



This final relation can also be derived directly without 
solving the differential equation at the surface: calculat- 
ing the value of the order parameter and its derivative at 
z = from (|49|) . one finds 



G 4 



Mo = 
dz 



2(1 



S 1 

G 4 



P) / dzf z n z , 
Jo 

f 00 



2(1 



r 



■G + (l-0) dzf zl j, z . 



The ratio fJ.' / fio then does not depend on the function 
Hz which has dropped out. As a result we recover 1)540. 
which also remains valid for the case of a bulk second- 
order phase transition (r = 0). In the limit B m 0, the 
boundary condition becomes /i = 0, since G ~ 0. Hence, 
the constant bulk solution goes through and the surface 
is not relevant, in agreement with the results of Sec. IV Al 
The analysis of the boundary value problem 145(1 with 
(|54|l follows the one in Ref. [l8|. Combining the expres- 
sion for the 'conserved energy' (originating from the bulk 
equation iPtSlO 



T 



(55) 



with the boundary condition 1(54(1 . we find an algebraic 
relation which determines the value of the order param- 
eter no at the surface, 



^G- 



T 



(56) 



The liquid surface fiQ = /xj = is always a solution and 
we deal with a continuous surface melting (O2 scenario) 
if it is the only one. Once a second solution with /io > is 
present, the surface undergoes a discontinuous (Oi) tran- 
sition, see Fig. Using the full expression for the po- 
tential on the RHS of 1(56(1. we find that 1(50(1 admits two 
solutions for large T (small fields B) and only the /io = 
solution for small temperatures (large magnetic fields), 
cf. the solid line in Fig. [3 Since both the continuous and 
the discontinuous melting scenarios arc present, a multi- 
critical point separating the two different kinds of transi- 
tions must exist. The equation which locates the critical 
point derives from ((56(1 by considering the quadratic ex- 
pansion of the potential Su!^ b (n)/T w 3(1 — c c )/i 2 (cf. 
Fig.EJ 



£{T mc ,B mc )G(B 

mc J 

^ (^rnc? ^mc) 



1 



(0.29 e d, 0.007 B x 



(57) 



For T < T mc the surface undergoes a continuous transi- 
tion and acts as a nucleus for the liquid phase, preventing 
the appearance of the solid metastable phase. On the 
other hand, for large temperatures T > T mc the order 
parameter at z = still undergoes a residual jump and 
the double-sided hysteretic behavior is restored, see Figs. 
Irtjlandim 
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FIG. 7: Residual value of the order parameter at the sur- 
face: analytic (solid line) from (I56H and full numerical solution 
(stars) see Section IVTI For T < T mc the surface undergoes a 
continuous transition (O2), whereas for T > T mc the order 
parameter at z = still exhibits a residual jump at melting. 



1. Multi- critical point at high-fields 



Zs 



FIG. 8: Sketch of a solid-liquid interface placed at a distance 
z a away from the surface. The surface destabilizing term 
E aul {(z s ) of Eq. 116011 favors the entry of the liquid-solid in- 
terface at the surface, producing a repulsive surface potential 
on the soliton. On the other hand, the energy cost to push 
the soliton into the system generates an attractive surface 
potential _E a i(z s ), cf. I6H . Depending on which term is dom- 
inant, either a continuous or a discontinuous surface melting 
transition is realized. 



At high magnetic fields the layers melt independently 
following a first-order type 2D-melting scenario. The or- 
der parameter /j,q in the topmost layer then undergoes 
a finite jump and the surface non- melting (Oi) scenario 
applies. The presence of a discontinuous 0% regime at 
high fields implies the existence of a second multi-critical 
point. While our analytical approach is not applica- 
ble anymore, since it is based on the gradient expan- 
sion which is not justified in this regime, numerically 
we have found clear indications of a finite jump at high 
fields {B « 10 B\). However, a more elaborate version 
of the DFT is required for an accurate determination of 
this multi-critical point. In particular, approaching the 
melting temperature of each 2D lattice, the higher-order 
peaks in the OCP correlation function c 2D (K) become 
important (K n > G). Hence, higher components in the 
Fourier expansion of the density have to be retained 40 , 
in order to obtain a more precise description of the high- 
field regime. 



C. Analysis of the continuous surface melting and 
of the multi-critical point 

The fact that a discontinuous bulk transition may turn 
continuous at the surface is somewhat non-trivial. To 
obtain a better insight into the mechanism of the con- 
tinuous surface melting, we present here an alternative 
analysis. We describe the surface-melting process in the 
language of the entry of the liquid through the boundary. 
The appearance of a liquid layer at the surface, while the 
bulk remains solid, implies the existence of an interface 
between the two phases, which is described by a soliton- 
like profile of the order parameter. The entry of the liquid 
can happen in two different ways: i) the soliton can slide 
smoothly from the boundary at T m (surface melting, O2) 
or ii) it starts entering the system but remains pinned 



near the surface at T m (surface non- melting, 0\). To 
distinguish these two scenarios we need to calculate the 
energy of the soliton as a function of its (half- height) po- 
sition z s away from the surface and the resulting pinning 
energy. We do this within the local theory described by 
the functional 

The saddle point equation of 1)58(1 reproduces the equa- 
tion of state l)45|l with the boundary condition (|54l) at 
z = (cf. Eqs. © and ©). 

As in the previous section, we concentrate on the case 
when the order parameter at the surface is small. Since 
the bulk solution remains valid also in the vicinity of the 
surface, we can study the problem within a variational 
scheme, taking the bulk soliton [i* (z a ) as a convenient 
variational function and displacing it rigidly at different 
distances z s from the surface ('si' stands for solid-liquid 
interface), see Pig. El 

We proceed with the estimate of the energy due to a 
solid- liquid interface at z s , which is given by 6u)[fig(z s )]. 
In the following, we want to calculate the soliton energy 
not only at the bulk transition temperature T m but also 
at temperatures T close by. However, for T ^ T m , the 
free energy 5u)[/j.f (z s )} is infinite due to the contribution 
from the bulk solid &J t 2 ° b (/x s ) ^ 0. In order to obtain 
a finite energy, which correctly estimates the soliton en- 
ergy E^ lf (z s ) in the presence of the surface we have to 
subtract this infinite contribution. We therefore define 

£T f (*s) = 8u[(4(z s )} ~ (L/d)6u;Z(»s,T) (59) 

where L — > 00 is thickness of the sample. Inserting the 
function fJ,f(z s ) in (|58|l . the third term is simple to cal- 
culate and yields the surface energy 

^HT^KfeM! (60) 
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The first two terms in i|58|l give the energy of the soliton 
in the semi- infinite system. For T < T m , the presence of 
a liquid layer of thickness ~ z s produces a contribution 
linear in z s , due to the difference between the liquid and 
solid free energies, — Sui^ b (fj, s ,T) > (we have written 
the additional argument T in 5L0^° b to make explicit its 
temperature dependence). Subtracting as in l|59|l the in- 
finite contribution from the solid phase, this energy be- 
comes — [(z s — z*)/d}5u>^ b (n s ,T), with z* an irrelevant 
constant. The remaining contribution to the soliton en- 
ergy is due to the interface E s i(z s ), i.e., the non-constant 
part of the order-parameter profile, and consists of both 
potential, &j s 2 ° b , and elastic, oc (/4) 2 , energies. Combin- 
ing these two terms, we find (see Eq. (BUI) ) 



E s i{z s ) 



IT <■» 
~~d 



dfi 



2^ 2 u D b (M) 



KM 

It 



T 



dfi 



E a i- y/6(l -Cc)T 



T 



(61) 



In the last line of l|61|) we have expanded the 2D potential 
for small values of fi. Since for a finite z s the soliton has 
not fully entered in the system, its energy is less than the 
total soliton energy E s \ in (|39|) due to the missing tail, see 
Pig. El One sees that, whereas -E s i penalizes the entrance 
of the soliton and thus favors the solid phase, E sur f pro- 
motes the formation of the liquid phase at the surface. 
Combining all terms, we obtain the soliton energy 1591) 



E^ rt (z s )^E sl + 
2d 



lT m -T 
P T m 



£G-^/6(l~c c ) {^(z s )} 2 , (62) 



where we have expanded the solid free energy term about 
T m , 6uZ(^,T) « d T 5Lu 2 ° b (fi S7 T m ){T-T m ) and we have 
used the definition of the latent heat I = T m As — 
-(p/d)T m d T 6u;™(p s ,T m ). 

Below melting, the term oc z s forbids the entry of the 
soliton into the bulk. At the melting transition T m , the 
third term plays the crucial role with a finite order pa- 
rameter at the surface //g (z s ) leading to a positive or a 
negative energy contribution depending on the sign of 
the prefactor oc [£G — ^/6(1 — c c )] . For a positive value 
(IG > 6(1 — c c )), the minimal energy is achieved by 
/Lig(z s ) = 0. As a consequence, the soliton enters com- 
pletely inside the sample at T m (surface melting). On 
the other hand, for £G < \/6(l — c c ), the entry of the 
interface costs a positive energy and, thus, the soliton 
remains pinned at the surface at T m . The two different 
behaviors are separated by the multicritical point, whose 
location is given by 



t{T mc , B mc )G(B mc ) = ^6(1 - So). (63) 
This equation coincides with l|57|) . 



The appearance of the multi-critical point (T mc , B mc ) 
along the melting line can be interpreted as a surface- 
depinning transition of the solid- liquid interface (soliton). 
To find the precise pinning location of the soliton in the 
0\ regime a more detailed analysis is required, account- 
ing for the higher order terms in (|61|l . We have calculated 
this potential numerically in both the surface melting and 
surface non-melting regimes, using the following scheme: 
we place a bulk soliton at different distances from the sur- 
face and evaluate numerically the total non-local energy 
(12911 as a function of z s . In Fig. we present two dif- 
ferent curves, one at T = 0.33 e^d > T mc and another at 
T = 0.28 £od < T mc . While for higher temperatures the 
potential exhibits a stable minimum close to the surface 
where the soliton remains pinned, upon decreasing the 
temperature the minimum moves deeper into the bulk 
and disappears altogether at T mc . 

When the surface melts continuously, we can charac- 
terize the transition by specific critical exponents, cf. Ref. 
ITsl In particular, we can look at how the soliton position 
z s depends on T. Here, we need the function /ig 1 (z s ), i.e., 
the relation between the value of the order parameter 
on the surface and the soliton position z s . For a soliton 
which is well inside the sample, /ig (z s ) can be approxi- 
mated as (cf. (0) 



(64) 



The position of the soliton is found from the minimum 
of (EIJ 



z s (T) w - 



1 



2G+v / 
0.5£|ln[(l 



:ln 



ld(T m - T) 



2^G+p£Tl(lG - v/6(l-c c )) 
T/T m )/(lG-tG\ m )]\. (65) 



Hence, the soliton slides into the system logarithmically 12 
with T — * T m . Next, we study how the residual order 
parameter on the surface goes to zero. Combining l|65|) 
together with <|64(l . we find that 



(66) 



These results are standard in the theory of surface melt- 
ing when only short-range interactions are present*^. 

For a continuous surface-melting transition the soliton 
propagates into the bulk at T m , leading, in a semi- infinite 
system, to the coexistence of the liquid and the solid. 
The interface deep in the bulk produces the maximum 
free energy cost, E^ rl (z s — + oo) = £? s i (see Sec. IIVB")l . 
This energy is only an interface energy and, hence, its 
contribution vanishes in the thermodynamic limit. In re- 
alistic finite systems, one has to account for the effect of 
the opposite surface. Approaching T m from below, both 
surfaces undergo a continuous melting transition. The 
two surfaces act like two nucleation points for the liquid 
phase, giving rise to two opposite solitons. The system 
is then composed of a sequence of liquid-solid-liquid re- 
gions. At T m the two solitons merge and the intermediate 
solid domain vanishes altogether. Hence, the solid cannot 
be overheated above the melting temperature. 
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FIG. 9: Energy as a function of the soliton position z a . Note 
the minimum for T = 0.33 eod > T mc (Oi, pinned soliton) 
which has disappeared at T = 0.28 eod < T mc (O2, depinned 
soliton) . 

VI. NUMERICAL ANALYSIS 

In order to check the accuracy of our analytical ap- 
proach, we have carried out a numerical solution of the 
saddle-point equations (|42[) . As a preliminary step, we 
write (|42[) in a slightly different form which is more con- 
venient for our numerical study. From 1|18|) . we obtain 
that at the saddle-point the solutions for £ z and [jl z are 
related by 

f°° dz' _ 

6 = y -^-Cz,z'fiz'- (67) 
Combining this expression with the relation (cf. (|24|l ) 

Mz = *'(&)/6, (68) 

where the function $ is defined in l)25|l. we obtain a set of 
integral equations which determine the order-parameter 
profile [i z 



/ d 2 Rg(R)exp 

J V 


- r°° dz' 

Jo ~ 


c z , z >[i z >g(R) 


d 2 K exp 

J V 


[J 


rco dz' 

3 -t c ~* 


x >(igig(R) 





where we have written $' explicitly. Our numerical solu- 
tion is based on the recursive solution of l|69|) . 

We first discretize the z axis in N = 1000 values 
{zi} with a fixed distance Zi — Zi—\ = Az = 0.04 A (for 
r > 40 we use a smaller step size Az = 0.01 A w d, 
since the soliton interface becomes sharper). We start 
from a constant solid phase and initialize the values of 
{fii} as [ii = fi Zi = [i s , for any i. Then, from Eq. 
1)6 7p. we derive the molecular field profile £j in corre- 
spondence to the values {zi}. We calculate the RHS 
of equation 1)68(1 for i = l,...,z max , while keeping the 
last values (i — i max + 1, . . . , 1000) unchanged, and ob- 
tain the new [if = 6$>'(£ Zi ). We compare {/i".} with 



{[i Zi } and if both the inequalities (1q — [1q < 10~ 5 and 
(1/N) Y^ii^z- ~ < 5 are satisfied, we accept {^".} 
as the order-parameter profile. Otherwise, the proce- 
dure is iterated recursively until a stable solution (fixed 
point) is reached. We take i max = 750; for this value 
we have checked that the connection between the numer- 
ical solution for i ^ « m ax and the constant bulk value 
for i > i max is smooth (we find a small jump at i max , 
(Mwx+i ~ ^bbI/Wb, ~ 10~ 5 )- Usually convergence is 
obtained after a reasonable number of iteration (< 100), 
however in the proximity of a continuous surface transi- 
tion the convergence becomes slow. This critical slowing 
down makes it difficult to track the sliding of the soliton 
inside the bulk. To avoid this problem, for these criti- 
cal cases we start our iteration from a more convenient 
initial state. Instead of initiating the profile in the bulk 
homogenous solid, we chose to start from a bulk soliton 
at the position which minimizes the total surface free en- 
ergy (cf. Fig. Eland discussion in the last section). In this 
case the convergence is extremely rapid. 

In Figs. Ell an d we show two different examples 
of the order-parameter profile for two different values of 
the temperature: i) T = 0.08 e d < T mc in Fig. fTUl 
where the surface undergoes a continuous transition and 
ii) T = 0.33 eod > T nlc in Fig. 1111 where a discontinu- 
ous surface transition takes place. For the O2 transition 
at T = 0.08 eod, we show the order-parameter profiles 
for different values of the magnetic field B, increasing 
from top to bottom. At low magnetic fields, the pro- 
file is almost constant (cf. the topmost line), since the 
surface kernel Soj s is negligible and the system is transla- 
tionally invariant. Going to larger values of the magnetic 
field, the modifications of the interlayer potential on the 
surface become relevant, reducing the value of the order 
parameter at the surface. Upon increasing the magnetic 
field further towards the thermodynamic melting transi- 
tion, the vortex density modulation becomes vanishingly 
small close to the surface. The numerical solution shows 
that the order parameter at the surface goes to zero (liq- 
uid) continuously. Hence, the surface assists the penetra- 
tion of the liquid phase by the formation of a quasi-liquid 
nucleus. This continuous transition on the surface elimi- 
nates the hysteresis above the melting transition, by pre- 
venting the appearance of a metastable overheated solid 
phase. 

At larger temperatures T = 0.33 eod > T mc (Fig. 
Illfl . the surface undergoes a discontinuous transition, al- 
though with a reduced jump in comparison with the bulk. 
Again starting from low magnetic fields (topmost line) 
the order parameter is constant, due to the smallness 
of the surface destabilizing potential. Increasing B, the 
value of [io decreases, still showing a larger suppression 
than in the bulk. However, the surface potential is not 
strong enough to push [io to zero and at melting (thick 
line) the surface still exhibits a finite order parameter. 
The transition to the homogeneous liquid phase then oc- 
curs via a finite global jump including the surface value 
[io as well. This discontinuous transition is compatible 
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FIG. 10: Numerical solutions of the order-parameter profile 
fjb z at the surface, for different values of the external magnetic 
field at T = 2e d/T = 25, corresponding to T = 0.08e d. 
At this temperature the bulk melts at Bm/Bx « 0.1528. 
The lines from top to bottom correspond to the magnetic 
fields with values B/Bx = 0.01, 0.05, 0.1, 0.125, 0.145, 0.15, 
B m /B\ - 10" 4 , Bm/B\ - 10" 5 , Bm/Bx - 10" 6 . While ap- 
proaching this value the soliton slides into the bulk (see the 
bottom three lines). The surface undergoes an O2 type tran- 
sition, i.e., no approaches zero (liquid phase) continuously. 



with the appearance of the metastable phase. Indeed, 
numerically it is possible to obtain a non-uniform solu- 
tion even at magnetic fields which are larger than the 
freezing field (lowest line in Fig. II If) . 

Finally, we check the accuracy of our analytical ap- 
proach in estimating the location of the multi-critical 
point (T mc , B mc ). We plot the residual value of the sur- 
face order parameter at melting as a function of temper- 
ature in Fig. [7| together with the solution of (|56l) . For a 
continuous surface melting transition the order parame- 
ter is exactly zero at melting. Numerically, we associate 
the value /io = to situations where the soliton poten- 
tial does not show a stable minimum as a function of z s 
(see Fig. Otherwise, we estimate the finite value of 
/io within the iterative solution of the saddle point equa- 
tions, which we have described above. The agreement of 
the numerical and analytical results is excellent, in par- 
ticular in the vicinity of the multi-critical point. 



VII. CONCLUSIONS 

In this paper, we have analyzed the impact of an ab- 
surface on the melting transition of the pancake vortex 
lattice. We have adapted the DFT-substrate approach 22 
to include the presence of the surface. We have found 
that for intermediate values of the magnetic field, the 
surface undergoes a continuous melting transition and 
assists a smooth propagation of the liquid into the bulk. 
This result reveals the origin of the asymmetric hystere- 
sis at the melting transition observed by Soibel et as 



12 zA 



FIG. 11: Numerical solutions of the order-parameter profile 
/i z near the surface at T = 0.33 eod (T = 6) for different values 
of the magnetic field. The freezing field is B^/Bx ~ 0.002396 
(thicker line). The lines from top to bottom correspond to 
the magnetic fields with values B/B x = 0, 0.0001, 0.0005, 
0.001, 0.0015, 0.002, 0.0022, 0.00235, £> m / Bx , B m /Bx + IQ~ 5 . 
The order parameter at the surface jumps discontinuously 
to zero; hence, the surface transition is in the class 0\. In 
this case, the surface does not preclude the appearance of 
the overheated solid. Numerically it is possible to obtain a 
non uniform solution even at magnetic fields B > B m larger 
than the freezing one (cf. lowest thin line), corresponding to 
a metastable configuration (overheated solid). 



the consequence of free surfaces which can act as nucle- 
ation sites for the liquid phase. Moreover, we have found 
that at low and high magnetic fields the surface transition 
turns discontinuous as in the bulk. The continuous and 
discontinuous transitions at the surface are separated by 
a multi-critical point. While the precise location of the 
high-field multi-critical point goes beyond the limits of 
validity of our analysis, we have determined the position 
of the low-field multi-critical point by means of an an- 
alytical solution of the DFT equations and numerically 
confirmed the result of our analysis. 

An effect which we have not included in our study is 
the impact of the reentrance of the melting line at low 
magnetic fields. At low magnetic fields the interaction 
between full vortex lines is strongly screened. This leads 
to a softening of the vortex lattice with exponentially 
small shear and compression moduli and to the reen- 
trance of the melting line towards small temperatures. 
The effect of a surface in this reentrant low-field regime 
is an interesting problem as well. The stray magnetic 
fields produce an algebraic interaction between the tips 
of the vortex finest instead of the exponentially small 
bulk interaction. This may lead to a highly unconven- 
tional scenario (proposed in Ref. l4l|) where the tips of 
the vortex lines on the surface are arranged in a regu- 
lar lattice while the vortex system is already melted in 
the bulk. Hence, in this case the surface plays a role 
which is opposite to the one played in the standard sur- 
face melting scenario, since it favors the appearance of 
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FIG. 12: Comparison of the soliton derived within the local 
approximation (thicker line) with the results of the full non- 
local theory for different values of Y = 6, 40, 80, 100, 130. For 
r = 6, upon rescaling z in units of £, the soliton collapses on 
the solution of the local approximation. Little deviations are 
present at F = 40, where the soliton is still well approximated 
by the gradient expansion. At higher F's the difference is ap- 
preciable: the full non-local theory leads to a shaper interface. 
Inset: plot of the ratio between £ and kernel extension 1/G+ 
as a function of F = 2eod/T. Assuming that the gradient ex- 
pansion is valid for £ C G+ ~ 4£G+ > 5, we obtain the limiting 
value F ~ 50, corresponding to T ~ 0.04 sod. 



where we have used I = \2a/dG\, cf. (EZJ - From the 
comparison of l|A3[) with l|A2|l it becomes clear that for 
small wave numbers k <C G+, i.e., for variations on scales 
larger than 1/G+, the gradient expansion approximates 
well the full non local theory. A convenient test is to 
compare the typical length of the soliton £ c 41 (cf. 
gTJ) with the size of the kernel 1/G+. If we choose 
£ C G+ w 5 as the limiting value for the validity of the 
expansion, we obtain that for T < 50 (corresponding to 
T > 0.04 eod or B < 0.5 B\) the gradient approximation 
is justified, see the inset of Fig. Expanding (|A3|) to 
second order, we find that the maximum error is only few 
percent, {k/G+f « 1/(£ G + ) 2 < 0.04. To confirm this 
result we have solved numerically the integral equation 
(|A1(I imposing the boundary conditions /Lt z _>_oo = and 
Mz— »oo = Ms- The results for different values of T are 
plotted in Fig. together with the soliton of the local 
theory (thick line). We have rescaled the z axis in units 
of the elastic length £. For T < 80 we obtain a good 
collapse of the data. On the other hand, for large T the 
kink in non-local theory is sharper than the soliton in 
the local theory and, therefore, the gradient expansion 
approximates poorly the exact result. 



the solid instead of the liquid. The impact of this 'sur- 
face solidification' on our results remains an interesting 
open question. 
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APPENDIX A: GRADIENT EXPANSION 



2. Linearised saddle-point equation 

For small values of /i z , we can expand the potential in 
the RHS of the equation IjAlfl , Hence, we write, cf. H46|l . 



1. Kernel in Fourier space 

In this Appendix we test the validity of the gradient 
expansion (|36|l for the DFT free energy in the low-field 
regime. We start from the expression of the bulk free- 
energy l|33[l and calculate the saddle point equation 



-6a 



dz' 



-/»-*'(/** - Av) = d^XM/T, (Al) 



which has to be compared with the corresponding equa- 
tion (|38|l originating from the approximated local theory. 
The difference between the two left-hand side terms of 
(|A1|) and (|38|l is more conveniently analyzed in Fourier 
space. In the case of the local theory, the second deriva- 
tive term is diagonal in Fourier space and its components 
read 



£ 2 k 2 ti k . 



(A2) 



In the full non-local theory, the left hand side of (|A1(I is 
also diagonal in Fourier space 



6a(/fc ~ //c=oW = 



1 + (k/G+) 



rMfc, 



(A3) 



M w «tW/ TRi6 ( 1 -^K 



(A4) 



where we consider the system at melting, i.e., c^° b = c c to 
connect with the discussion in Sec. IV Bl Inserting (|A4|I 
and l|A3(l in I|A1(I . we obtain the equation of 'motion' in 
Fourier space 



-ek 1 



6(1 - Cc) 



(£G + y 



fi k = 6(1 - c c )fi k . 



(A5) 



We can transform this equation back to real space 



(1 + r)£ 2 fi" = 6(1 - c c )fi z 



(A6) 



where we have defined the parameter r = dG + (l — 
c c )/2a = 6(1 - c c )/(G+£) 2 , as in Eq. (JSTJ. The higher 
derivatives which are neglected in the gradient expan- 
sion produce a small renormalization of the elastic term 
(1+r) « 1. 
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